/*
* INSANE - Interactive Structural Analysis Environment
*
* Copyright (C) 2003-2005
* Universidade Federal de Minas Gerais
* Escola de Engenharia
* Departamento de Engenharia de Estruturas
* 
* Author's email :    insane@dees.ufmg.br
* Author's website :  http://www.dees.ufmg.br/insane
* 
* This program is free software; you can redistribute it and/or
* modify it under the terms of the GNU General Public License
* as published by the Free Software Foundation; either version 2
* of the License, or any later version.
* 
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
* GNU General Public License for more details.
* 
* You should have received a copy of the GNU General Public License
* along with this program; if not, write to the Free Software
* Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
*/
package br.ufmg.dees.insane.model.femModel.problemDriver;

import java.io.Serializable;

import br.ufmg.dees.insane.model.femModel.element.Element;
import br.ufmg.dees.insane.model.femModel.element.ThinPlateElement;
import br.ufmg.dees.insane.util.IMatrix;

/** 
 * A class representing CowperProblemDriver
 *
 * @author Samir
 * @since 15/09/2006
 */
public class CowperProblemDriver extends KirchhoffThinPlate implements Serializable{

	private static final long serialVersionUID = 1L;
	double a = 0.0;
	double b = 0.0;
	double c = 0.0;

	 /* (non-Javadoc)
	 * @see br.ufmg.dees.insane.model.femModel.problemDriver.ProblemDriver#getC(br.ufmg.dees.insane.model.femModel.element.Element)
	 */
	public IMatrix getC(Element e) {
	    int numberOfDegreesOfFreedom = ((ThinPlateElement)e).getNumberOfDegreesOfFreedom();
		IMatrix finalStiffnessMatrix = new IMatrix(numberOfDegreesOfFreedom,numberOfDegreesOfFreedom);
	    finalStiffnessMatrix.setZero();
	    IMatrix cartesianNodalCoordsMatrix = e.getCartesianNodalCoordsMatrix();
	    
	    double x1 = cartesianNodalCoordsMatrix.getElement(0,0);
	    double x2 = cartesianNodalCoordsMatrix.getElement(1,0);
	    double x3 = cartesianNodalCoordsMatrix.getElement(2,0);

	    double y1 = cartesianNodalCoordsMatrix.getElement(0,1);
	    double y2 = cartesianNodalCoordsMatrix.getElement(1,1);
	    double y3 = cartesianNodalCoordsMatrix.getElement(2,1);
	    
		double r = Math.sqrt(Math.pow((x2-x1),2)+Math.pow((y2-y1),2));
		
		this.a = ((x2-x3)*(x2-x1)+(y2-y3)*(y2-y1))/r;
		this.b = ((x3-x1)*(x2-x1)+(y3-y1)*(y2-y1))/r;
		this.c = ((x2-x1)*(y3-y1)-(x3-x1)*(y2-y1))/r;
	    	    
	    double t = ((ThinPlateElement)e).getThickness();
	    double E = ((ThinPlateElement)e).getYoungsModulus();
	    double v = ((ThinPlateElement)e).getPoissonsRatio();
	    
	    double D = (E*t*t*t)/(12*(1-v*v));
	    
	    IMatrix T = this.T();
	    IMatrix R = this.R(x1,x2,y1,y2);
	    IMatrix k = this.k(v);
	    
	    IMatrix T2R = new IMatrix(20,18);
	    T2R.setZero();
	    T2R.mul(this.T2(T),R);
	    	    
	    IMatrix K1 = new IMatrix(T2R.getNumCol(),k.getNumCol());
	    K1.mulTransposeLeft(T2R,k);
	    
	    finalStiffnessMatrix.mul(K1,T2R);
	    finalStiffnessMatrix.setScale(D);
	    
	    return(finalStiffnessMatrix);
	}
	
	private int[] m(){
		int[] m = new int[20];
		m[0] = 0;	m[1] = 1;	m[2] = 0;	m[3] = 2;
		m[4] = 1;	m[5] = 0;	m[6] = 3;	m[7] = 2;
		m[8] = 1;	m[9] = 0;	m[10] = 4;	m[11] = 3;
		m[12] = 2;	m[13] = 1;	m[14] = 0;	m[15] = 5;
		m[16] = 3;	m[17] = 2;	m[18] = 1;	m[19] = 0;
		
		return m;
	}

//	private IVector m(){
//		IVector m = new IVector(20);
//		m.setElement(0,0);	m.setElement(1,1);	m.setElement(2,0);	m.setElement(3,2);
//		m.setElement(4,1);	m.setElement(5,0);	m.setElement(6,3);	m.setElement(7,2);
//		m.setElement(8,1);	m.setElement(9,3);	m.setElement(10,4);	m.setElement(11,3);
//		m.setElement(12,2);	m.setElement(13,1);	m.setElement(14,0);	m.setElement(15,5);
//		m.setElement(16,3);	m.setElement(17,2);	m.setElement(18,1);	m.setElement(19,0);
//		
//		return m;
//	}
	
	private int[] n(){
		int[] n = new int[20];
		n[0] = 0;	n[1] = 0;	n[2] = 1;	n[3] = 0;
		n[4] = 1;	n[5] = 2;	n[6] = 0;	n[7] = 1;
		n[8] = 2;	n[9] = 3;	n[10] = 0;	n[11] = 1;
		n[12] = 2;	n[13] = 3;	n[14] = 4;	n[15] = 0;
		n[16] = 2;	n[17] = 3;	n[18] = 4;	n[19] = 5;
		
		return n;
	}
	
//	private IVector n(){
//		IVector n = new IVector(20);
//		n.setElement(0,0);	n.setElement(1,0);	n.setElement(2,1);	n.setElement(3,0);
//		n.setElement(4,1);	n.setElement(5,2);	n.setElement(6,0);	n.setElement(7,1);
//		n.setElement(8,2);	n.setElement(9,0);	n.setElement(10,0);	n.setElement(11,1);
//		n.setElement(12,2);	n.setElement(13,3);	n.setElement(14,4);	n.setElement(15,0);
//		n.setElement(16,2);	n.setElement(17,3);	n.setElement(18,4);	n.setElement(19,5);
//		
//		return n;
//	}
	
	private double F(int m, int n){
		double F = 0.0;
//		double f1 = 0.0;  double f21 = 0.0;  double f22 = 0.0;  double f3 = 0.0;  double f4 = 0.0;
//		int n1 = n + 1;
//		int m1 = m + 1;
//		
//		f1 = Math.pow(this.c, n1);
//		f21 = Math.pow(this.a, m1);
//		f22 = Math.pow(-this.b, m1);
//		f3 = this.factorial(m)*this.factorial(n);
//		f4 = this.factorial((m+n+2));
//		
//		F = f1*(f21-f22)*f3/f4;
		
		F = (Math.pow(this.c,(n+1))*(Math.pow(this.a,(m+1))-(Math.pow((-this.b),(m+1))))*this.factorial(m)*this.factorial(n))/this.factorial((m+n+2));
				
		return F;
	}
	
	private IMatrix k(double v){
		IMatrix k = new IMatrix(20,20);
		k.setZero();
		int[] m = this.m();
		int[] n = this.n();
		
		for(int i = 0; i < 20; i++){
			for(int j = 0; j < 20; j++){
				double kij = 0.0;  double k1 = 0.0;  double k2 = 0.0;  double k3 = 0.0;   double k4 = 0.0;	
				int mi = m[i];  int mj = m[j];  int ni = n[i];	int nj = n[j];
				
				if(mi<2 || mj<2){
					k1 = 0;
				} else {
					k1 = mi*mj*(mi-1)*(mj-1)*this.F(mi+mj-4, ni+nj);
				}
				if(ni<2 || nj<2){
					k2 = 0;
				} else {
					k2 = ni*nj*(ni-1)*(nj-1)*this.F(mi+mj, ni+nj-4);
				}
				k3 = 2*(1-v)*mi*mj*ni*nj + v*mi*nj*(mi-1)*(nj-1)+v*mj*ni*(mj-1)*(ni-1);
				if(k3==0){
					k4 = 0;
				} else {
					k4 = this.F(mi+mj-2, ni+nj-2);
				}
				kij = k1 + k2 + (k3*k4); 

//				kij = m.getElement(i)*m.getElement(j)*(m.getElement(i)-1)*(m.getElement(j)-1)*
//					  this.F((int)(m.getElement(i)+m.getElement(j)-4),(int)(n.getElement(i)+n.getElement(j)))+
//					  n.getElement(i)*n.getElement(j)*(n.getElement(i)-1)*(n.getElement(j)-1)*
//					  this.F((int)(m.getElement(i)+m.getElement(j)),(int)(n.getElement(i)+n.getElement(j)-4))+
//					  (2*(1-v)*m.getElement(i)*m.getElement(j)*n.getElement(i)*n.getElement(j)+
//					  v*m.getElement(i)*n.getElement(j)*(m.getElement(i)-1)*(n.getElement(j)-1)+
//					  v*m.getElement(j)*n.getElement(i)*(m.getElement(j)-1)*(n.getElement(i)-1))*
//					  this.F((int)(m.getElement(i)+m.getElement(j)-2),(int)(n.getElement(i)+n.getElement(j)-2));
				k.setElement(i,j,kij);
			}
		}
			
		return k;
	}
	
	private IMatrix T(){
		IMatrix T = new IMatrix(20,20);
		T.setZero();
		
		T.setElement(0,0,1);  T.setElement(0,1,-this.b);  T.setElement(0,3,Math.pow(this.b,2));
		T.setElement(0,6,-Math.pow(this.b,3));  T.setElement(0,10,Math.pow(this.b,4));  T.setElement(0,15,-Math.pow(this.b,5));
		
		T.setElement(1,1,1);  T.setElement(1,3,-2*this.b);  T.setElement(1,6,3*Math.pow(this.b,2));
		T.setElement(1,10,-4*Math.pow(this.b,3));  T.setElement(1,15,5*Math.pow(this.b,4));
		
		T.setElement(2,2,1);  T.setElement(2,4,-this.b);  T.setElement(2,7,Math.pow(this.b,2));
		T.setElement(2,11,-Math.pow(this.b,3));
		
		T.setElement(3,3,2);  T.setElement(3,6,-6*this.b);  T.setElement(3,10,12*Math.pow(this.b,2));
		T.setElement(3,15,-20*Math.pow(this.b,3));
		
		T.setElement(4,4,1);  T.setElement(4,7,-2*this.b);  T.setElement(4,11,3*Math.pow(this.b,2));
		
		T.setElement(5,5,2);  T.setElement(5,8,-2*this.b);  T.setElement(5,12,2*Math.pow(this.b,2));
		T.setElement(5,16,-2*Math.pow(this.b,3));
		
		T.setElement(6,0,1);  T.setElement(6,1,a);  T.setElement(6,3,Math.pow(this.a,2));
		T.setElement(6,6,Math.pow(this.a,3));  T.setElement(6,10,Math.pow(this.a,4));  T.setElement(6,15,Math.pow(this.a,5));
		
		T.setElement(7,1,1);  T.setElement(7,3,2*this.a);  T.setElement(7,6,3*Math.pow(this.a,2));
		T.setElement(7,10,4*Math.pow(this.a,3));  T.setElement(7,15,5*Math.pow(this.a,4));
		
		T.setElement(8,2,1);  T.setElement(8,4,this.a);  T.setElement(8,7,Math.pow(this.a,2));
		T.setElement(8,11,Math.pow(this.a,3));
		
		T.setElement(9,3,2);  T.setElement(9,6,6*this.a);  T.setElement(9,10,12*Math.pow(this.a,2));
		T.setElement(9,15,20*Math.pow(this.a,3));
		
		T.setElement(10,4,1);  T.setElement(10,7,2*this.a);  T.setElement(10,11,3*Math.pow(this.a,2));
		
		T.setElement(11,5,2);  T.setElement(11,8,2*this.a);  T.setElement(11,12,2*Math.pow(this.a,2));
		T.setElement(11,16,2*Math.pow(this.a,3));
		
		T.setElement(12,0,1);  T.setElement(12,2,this.c);  T.setElement(12,5,Math.pow(this.c,2));
		T.setElement(12,9,Math.pow(this.c,3));  T.setElement(12,14,Math.pow(this.c,4));  T.setElement(12,19,Math.pow(this.c,5));
		
		T.setElement(13,1,1);  T.setElement(13,4,this.c);  T.setElement(13,8,Math.pow(this.c,2));
		T.setElement(13,13,Math.pow(this.c,3));  T.setElement(13,18,Math.pow(this.c,4));
		
		T.setElement(14,2,1);  T.setElement(14,5,2*this.c);  T.setElement(14,9,3*Math.pow(this.c,2));
		T.setElement(14,14,4*Math.pow(this.c,3));  T.setElement(14,19,5*Math.pow(this.c,4));
		
		T.setElement(15,3,2);  T.setElement(15,7,2*this.c);  T.setElement(15,12,2*Math.pow(this.c,2));
		T.setElement(15,17,2*Math.pow(this.c,3));
		
		T.setElement(16,4,1);  T.setElement(16,8,2*this.c);  T.setElement(16,13,3*Math.pow(this.c,2));
		T.setElement(16,18,4*Math.pow(this.c,3));
		
		T.setElement(17,5,2);  T.setElement(17,9,6*this.c);  T.setElement(17,14,12*Math.pow(this.c,2));
		T.setElement(17,19,20*Math.pow(this.c,3));
		
		T.setElement(18,15,5*Math.pow(this.a,4)*this.c);  
		T.setElement(18,16,(3*Math.pow(this.a,2)*Math.pow(this.c,3)-2*Math.pow(this.a,4)*this.c));
		T.setElement(18,17,(-2*this.a*Math.pow(this.c,4)+3*Math.pow(this.a,3)*Math.pow(this.c,2)));
		T.setElement(18,18,(Math.pow(this.c,5)-4*Math.pow(this.a,2)*Math.pow(this.c,3)));
		T.setElement(18,19,5*this.a*Math.pow(this.c,4));
		
		T.setElement(19,15,5*Math.pow(this.b,4)*this.c);  
		T.setElement(19,16,(3*Math.pow(this.b,2)*Math.pow(this.c,3)-2*Math.pow(this.b,4)*this.c));
		T.setElement(19,17,(2*this.b*Math.pow(this.c,4)-3*Math.pow(this.b,3)*Math.pow(this.c,2)));
		T.setElement(19,18,(Math.pow(this.c,5)-4*Math.pow(this.b,2)*Math.pow(this.c,3)));
		T.setElement(19,19,-5*this.b*Math.pow(this.c,4));

		return T;
	}
	
	private IMatrix R(double x1, double x2, double y1, double y2){
		IMatrix R = new IMatrix(18,18);
		R.setZero();
		
		double r = Math.sqrt(Math.pow((x2-x1),2)+Math.pow((y2-y1),2));
		double cos = (x2-x1)/r;
		double sen = (y2-y1)/r;
		
		R.setElement(0,0,1);
		R.setElement(1,1,cos);  R.setElement(1,2,sen);
		R.setElement(2,1,-sen);  R.setElement(2,2,cos);
		R.setElement(3,3,Math.pow(cos,2));  R.setElement(3,4,2*sen*cos);  R.setElement(3,5,Math.pow(sen,2));
		R.setElement(4,3,-sen*cos);  R.setElement(4,4,(Math.pow(cos,2)-Math.pow(sen,2)));  R.setElement(4,5,sen*cos);
		R.setElement(5,3,Math.pow(sen,2));  R.setElement(5,4,-2*sen*cos);  R.setElement(5,5,Math.pow(cos,2));
		
		R.setElement(6,6,1);
		R.setElement(7,7,cos);  R.setElement(7,8,sen);
		R.setElement(8,7,-sen);  R.setElement(8,8,cos);
		R.setElement(9,9,Math.pow(cos,2));  R.setElement(9,10,2*sen*cos);  R.setElement(9,11,Math.pow(sen,2));
		R.setElement(10,9,-sen*cos);  R.setElement(10,10,(Math.pow(cos,2)-Math.pow(sen,2)));  R.setElement(10,11,sen*cos);
		R.setElement(11,9,Math.pow(sen,2));  R.setElement(11,10,-2*sen*cos);  R.setElement(11,11,Math.pow(cos,2));
		
		R.setElement(12,12,1);
		R.setElement(13,13,cos);  R.setElement(13,14,sen);
		R.setElement(14,13,-sen);  R.setElement(14,14,cos);
		R.setElement(15,15,Math.pow(cos,2));  R.setElement(15,16,2*sen*cos);  R.setElement(15,17,Math.pow(sen,2));
		R.setElement(16,15,-sen*cos);  R.setElement(16,16,(Math.pow(cos,2)-Math.pow(sen,2)));  R.setElement(16,17,sen*cos);
		R.setElement(17,15,Math.pow(sen,2));  R.setElement(17,16,-2*sen*cos);  R.setElement(17,17,Math.pow(cos,2));

		return R;
	}
	
	private IMatrix T2(IMatrix T){
		IMatrix T2 = new IMatrix(20,18);
		T2.setZero();
		
		IMatrix invT = T.pseudoInverse();
		
		for(int i = 0; i < 20; i++){
			for(int j = 0; j < 18; j++){
				T2.setElement(i,j,invT.getElement(i,j));
			}
		}
		
		return T2;
	}

	private int factorial(int num){
		int f = 0;
		int cont = num - 1;
		f = num;
		for(int i = 0; i < cont; i++){
			num = num - 1;
			f = f*num;
		}
		if(num==0 || num < 0){
			f = 1;
		}
		
		return f;
	}

	public String getUnknownDualVariableLabel() {
		// TODO Auto-generated method stub
		return null;
	}
	
//	private double pow(double a, int b){
//		double r = 0.0;
//		r = a;
//		if(b > 0){
//			for (int i = 1; i < b; i++) {
//				r = r*r;
//			}
//		} else if(b < 0){
//			for (int i = 1; i < -b; i++) {
//				r = r*r;
//			}
//			r = 1.0/r;
//		} else if(b == 0){
//			r = 1.0;
//		}
//		return r;
//	}
}
